program SS_JPEM_Main

!====================================================================================================
!
! "Borrowing Premia in Unsecured Credit Markets"
! by Kyle Dempsey and Felicia Ionescu
! 01.24.2024
!
! This code solves the extended version of the model constructed for revision at the JPE:Macro.
! The main additions relative to the original submission are:
! 	- long term debt
! 	- stochastic choice on default AND borrowing / saving actions 
! 	- richer default specifications (liquidation, reorganization, payment plan)
!
! Code by: Kyle Dempsey
!
! NOTES:
! - all parameters, program objects, read-in files, etc. defined at top of module file.
! - original version of code would have "SingleUpdateMethod" set to 0
!
!===================================================================================================

use omp_lib
use kind_mod
use misc_utilities
use SS_JPEM_Module

implicit none

real(rk) :: tmp, afval, aw, prob, lpval, CriterionValue, KBar	 ! average / constant borrowing premium for fixed premium case
integer(ik) :: ai, e1i, e2i, e3i, i, afi, tmpi



select case (ModelVersion)
	case(0)
		print*, 'BASELINE'
	case(1)	
		print*, 'SHORT TERM DEBT MODEL'
		CaseDir = ShortTermDir 
	case(2)
		print*, 'LOW INTERMEDIATION COST'
		CaseDir = LowIntCostDir 
	case(3)
		print*, 'FIXED RECOVERY: CH 7 ONLY'
		CaseDir = FixedRecoveryDir 
	case(4)
		print*, 'BANKRUPTCY ONLY'
		CaseDir = BKOnlyDir 
end select 

InDir = trim(MasterInDir) // trim(CaseDir)
OutDir = trim(MasterOutDir) // trim(CaseDir)
PanelDir = trim(MasterOutDir) // trim(CaseDir) // trim(MasterPanelDir)
print*, ' '



!-------------------------------!
! PARAMETERS AND TARGET MOMENTS !
!-------------------------------!
! EXTERNAL PARAMETERS 
812 format(a30, a10, 3f8.3)
print*, ' '
print*, ' EXTERNAL PARAMETERS'
print*, ' '
print*, 'technology and preferences'
print 812, 'capital share','alpha', alpha
print 812, 'minimum consumption','cmin', cmin
print 812, 'risk aversion','crra', crra
print*, ' '
print*, 'asset market'
print 812, 'Macaulay Duration', 'D', Macaulay
print 812, 'risk-free price normalization', 'qBar', qBar
print 812, 'fixed cost of lending', 'F', LoanFixedCost
print*, ' '
print*, 'legal'
print 812, 'exclusion, Ch 7','1/theta_7', 1.0_rk / theta_7
print 812, 'exclusion, Ch 13','1/theta_13', 1.0_rk / theta_13
print 812, 'discharge rate, Ch 13','xi_13', xi_13
print 812, 'max garnishment, Ch 13','lambdaBar', lambdaBar
print 812, 'filing cost', 'psi', psi
print*, ' '


allocate(ExternalParameters(15))


allocate(InternalParameters(ncal))
InternalParameters(1) = beta
InternalParameters(2) = stigma_7
InternalParameters(3) = stigma_13
InternalParameters(4) = stigma_R
InternalParameters(5) = ev_nu
InternalParameters(6) = ev_zeta
InternalParameters(7) = ev_eta
InternalParameters(8) = IntermediationCost
InternalParameters(9) = LoanFixedCost
InternalParameters(10) = alpha_R


allocate(InternalParameters_LB(ncal),InternalParameters_UB(ncal))
InternalParameters_LB(1) = beta_LB; InternalParameters_UB(1) = beta_UB
InternalParameters_LB(2) = stigma_7_LB; InternalParameters_UB(2) = stigma_7_UB
InternalParameters_LB(3) = stigma_13_LB; InternalParameters_UB(3) = stigma_13_UB
InternalParameters_LB(4) = stigma_R_LB; InternalParameters_UB(4) = stigma_R_UB
InternalParameters_LB(5) = ev_nu_LB; InternalParameters_UB(5) = ev_nu_UB
InternalParameters_LB(6) = ev_zeta_LB; InternalParameters_UB(6) = ev_zeta_UB
InternalParameters_LB(7) = ev_eta_LB; InternalParameters_UB(7) = ev_eta_UB
InternalParameters_LB(8) = IntermediationCost_LB; InternalParameters_UB(8) = IntermediationCost_UB
InternalParameters_LB(9) = LoanFixedCost_LB; InternalParameters_UB(9) = LoanFixedCost_UB
InternalParameters_LB(10) = alpha_R_LB; InternalParameters_UB(10) = alpha_R_UB


! EMPIRICAL TARGETS
allocate(Data_TargetMoments(ntarg))
print*, ' '
print*, 'TARGETING DATA MOMENTS'
Data_TargetMoments(1) = Data_BankruptcyRate
Data_TargetMoments(2) = Data_Ch7Share
Data_TargetMoments(3) = Data_RecoveryRate
Data_TargetMoments(4) = Data_ChargeoffRate
Data_TargetMoments(5) = Data_FractionInDebt
Data_TargetMoments(6) = Data_DebtToIncomeRatio
Data_TargetMoments(7) = Data_AverageLoanSpread
Data_TargetMoments(8) = Data_SpreadIntercept
Data_TargetMoments(9) = Data_SpreadSlope
print*, ' '

allocate(TargetWeights(ntarg))
TargetWeights(1) = Weight_BankruptcyRate
TargetWeights(2) = Weight_Ch7Share
TargetWeights(3) = Weight_RecoveryRate
TargetWeights(4) = Weight_ChargeoffRate
TargetWeights(5) = Weight_FractionInDebt
TargetWeights(6) = Weight_DebtToIncomeRatio
TargetWeights(7) = Weight_AverageLoanSpread
TargetWeights(8) = Weight_SpreadIntercept
TargetWeights(9) = Weight_SpreadSlope

211 format(a40, 2a8)
212 format(a40, 2f8.2)
print*, ' '
print*, ' '
print*, ' '
print*, ' '
print*, ' '
print 211, ' target and model moments (all in %)   ', 'data', 'weight'
print 212, ' bankruptcy rate                       ', Data_TargetMoments(1), TargetWeights(1)
print 212, ' Ch 7 share of bankruptcy              ', Data_TargetMoments(2), TargetWeights(2)
print 212, ' recovery rate                         ', Data_TargetMoments(3), TargetWeights(3)
print 212, ' chargeoff rate                        ', Data_TargetMoments(4), TargetWeights(4)
print 212, ' fraction in debt                      ', Data_TargetMoments(5), TargetWeights(5)
print 212, ' debt to income ratio                  ', Data_TargetMoments(6), TargetWeights(6)
print 212, ' average interest rate spread on loans ', Data_TargetMoments(7), TargetWeights(7)
print 212, ' spread profile intercept              ', Data_TargetMoments(8), TargetWeights(8)
print 212, ' spread profile slope                  ', Data_TargetMoments(9), TargetWeights(9)

!-------!
! GRIDS !
!-------!
print*, ' '
print*, 'creating grids...'
! ASSETS 
!(coarse, negative and positive regions both log spaced, then combined)
allocate(aGrid(na), aGridNeg(naNeg), aGridPos(naPos))
call logspace(-amax_neg - 1.0_rk * (-amax_neg) + 1.0_rk, -amin_neg - 1.0_rk * (-amax_neg) + 1.0_rk, naNeg, aGridNeg)
aGridNeg = aGridNeg - amax_neg - 1.0_rk  ! comes out positive and directionally wrong!
do ai = 1, naNeg
    aGrid(ai) = - aGridNeg(naNeg + 1 - ai)
end do
aGridNeg = aGrid(1:naNeg)
call logspace( amin_pos - 1.0_rk * amin_pos + 1.0_rk, amax_pos - 1.0_rk * amin_pos + 1.0_rk, naPos, aGridPos)	
aGridPos = aGridPos - 1.0_rk
aGrid(aZeroInd:na) = aGridPos


! EARNINGS
! NB: combination of e2 / e3 / beta into a single combined state must be done in the model solution step 
! for a given set of parameters, because these parameters include beta levels and transitions 
432 format(i5,20f10.4)

! PERMANENT COMPONENT
allocate(e1Grid(ne1))
call linspace(-nsd1 * sigma1, nsd1 * sigma1, ne1, e1Grid)	! log(e1) grid 
allocate(e1Share(ne1)) 
e1Share = normal_pdf(e1Grid, 0.0_rk, sigma1)				! e1 prob 
e1Grid = exp(e1Grid)										! convert to levels from logs 
e1Share = e1Share / sum(e1Share)							! normalize from continuous PDF analog
print*, 'PERMANENT EARNINGS'
do e1i = 1,ne1
	print 432, e1i, e1Grid(e1i), e1Share(e1i)
end do 
print*, ' '


! PERMANENT COMPONENT 
allocate(e2Grid_Orig(ne2), e2Trans(ne2, ne2))
call tauchen(rho2, 0.0_rk, sigma2, ne2, nsd2, e2Grid_Orig, e2Trans)
e2Grid_Orig = exp(e2Grid_Orig)								! convert to levels from logs 
allocate(e2Stat(ne2))
e2Stat = e2Trans(1,:)
do i = 1,100000
    e2Stat = matmul(e2Stat, e2Trans)
end do
print*, 'PERSISTENT EARNINGS'
do e2i = 1,ne2
	print 432, e2i, e2Grid_Orig(e2i), e2Stat(e2i), e2Trans(e2i,:)
end do 
print*, ' '


! TRANSITORY COMPONENT
allocate(e3Grid_Orig(ne3))
call linspace(-nsd3 * sigma3, nsd3 * sigma3, ne3, e3Grid_Orig)	! log(e3) grid 
allocate(e3Prob(ne3)) 
e3Prob = normal_pdf(e3Grid_Orig, 0.0_rk, sigma3)				! e3 prob 
e3Grid_Orig = exp(e3Grid_Orig)									! convert to levels from logs 
e3Prob = e3Prob / sum(e3Prob)									! normalize from continuous PDF analog
print*, 'TRANSITORY EARNINGS'
do e3i = 1,ne3
	print 432, e3i, e3Grid_Orig(e3i), e3Prob(e3i)
end do 
print*, ' '


! create cumulative versions for e2 and e3 
allocate(cum_e2Trans(ne2, ne2), cum_e3Prob(ne3))
call CumulativeTransitionMatrix(cum_e2Trans, e2Trans, ne2)
call CumulativeProbabilityVector(cum_e3Prob, e3Prob, ne3)


! compute aggregate labor productivity
allocate(lpMat(ne1, ne2, ne3))
AggregateLabor = 0.0_rk
do e3i = 1, ne3
do e2i = 1, ne2
do e1i = 1, ne1
	lpVal = e1Grid(e1i) * e2Grid_Orig(e2i) * e3Grid_Orig(e3i)
	lpMat(e1i, e2i, e3i) = lpVal
    AggregateLabor = AggregateLabor + e1Share(e1i) * e2Stat(e2i) * e3Prob(e3i) * lpVal
end do
end do
end do

! compute prices and things implied by the labor process given the K-Y ratio 
KBar = AggregateLabor * (Data_KYRatio ** (1.0 / (1.0_rk - alpha)))
wage = (1.0_rk - alpha) * (KBar / AggregateLabor) ** alpha 

! garnishment
allocate(lambdaGrid(nl))
call linspace(0.0_rk, lambdaBar, nl, lambdaGrid)

!----------------------------------------!
! INITIAL GUESSES OF EQUILIBRIUM OBJECTS !
!----------------------------------------!
print*, 'initializing equilibrium objects...'
print*, 'value functions...'
! ex ante value functions
allocate(W0(na, ni, ne1), WN0(na, ni, ne1), W7(naPos, ni, ne1), W13(naPos, ni, ne1, nl))
if (WReadIn .eq. 1_ik) then 
	call reading3d(W0, na, ni, ne1, 'W0', InDir)
	call reading3d(WN0, na, ni, ne1, 'WN0', InDir)
	call reading3d(W7, naPos, ni, ne1, 'W7', InDir)
	call reading4d(W13, naPos, ni, ne1, nl, 'W13', InDir)
else
	W0 = 0.0_rk 	! no EGM in this version so this is a fine guess.
	W7 = 0.0_rk
	W13 = 0.0_rk 
	WN0 = 0.0_rk 
end if 

! a' decisions
print*, 'saving decisions...'
allocate(pA0(na,na,ni,ne1), pA7(naPos, naPos, ni, ne1), pA13(naPos, naPos, ni, ne1, nl))
if (pAReadIn .eq. 1_ik) then 
	call reading4d(pA0, na, na, ni, ne1, 'pA0', InDir)
	call reading4d(pA7, naPos, naPos, ni, ne1, 'pA7', InDir)
	call reading5d(pA13, naPos, naPos, ni, ne1, nl, 'pA13', InDir)
else
	pA0 = 0.0_rk
	do ai = 1,na
		pA0(aZeroInd, ai, :, :) = 1.0_rk 	! just assume everyone goes to zero
		if (ai .ge. aZeroInd) then 
			pA7(1, aiPositive(ai), :, :) = 1.0_rk
			pA13(1, aiPositive(ai), :, :,:) = 1.0_rk
		end if 
	end do 
end if 

! default and default type decisions
print*, 'default decisions...'
allocate(pD(naNeg,ni,ne1), pD7(naNeg, ni, ne1), pD13(naNeg, ni, ne1), pDR(naNeg, ni, ne1)) 
if (pDReadIn .eq. 1_ik) then 
	call reading3d(pD, naNeg, ni, ne1, 'pD', InDir)
	call reading3d(pD7, naNeg, ni, ne1, 'pD7', InDir)
	call reading3d(pD13, naNeg, ni, ne1, 'pD13', InDir)
	call reading3d(pDR, naNeg, ni, ne1, 'pDR', InDir)
else
	pD = 0.0_rk 	! just assume no one defaults 
	pD7 = 1.0_rk 	! just assume everyone liquidates conditional on default 
	pD13 = 0.0_rk
	pDR = 0.0_rk
end if 

! stationary distribution
print*, 'distributions...'
allocate(mu0(na,ni,ne1), mu7(naPos, ni, ne1), mu13(naPos, ni, ne1, nl))
if (muReadIn .eq. 1_ik) then 
	call reading3d(mu0, na, ni, ne1, 'mu0', InDir)
	call reading3d(mu7, naPos, ni, ne1, 'mu7', InDir)
	call reading4d(mu13, naPos, ni, ne1, nl, 'mu13', InDir)
end if 


! loan prices
print*, 'loan prices...'
allocate(q(naNeg,na,ni,ne1))
allocate(XD(naNeg,ni,ne1), XN(naNeg,ni,ne1))
if (qReadIn .eq. 1_ik) then 
	call reading4d(q, naNeg, na, ni, ne1, 'q', InDir)
	call reading3d(XD, naNeg, ni, ne1, 'XD', InDir)
	call reading3d(XN, naNeg, ni, ne1, 'XN', InDir)
else
	q = qBar	! set to risk-free price
	XD = 0.0_rk 
	XN = 1.0_rk 
end if 
print*, ' '

!==========================!
!                          !
! SOLVE OR CALIBRATE MODEL !
!                          !
!==========================!
print*, 'completing allocations...'
! idiosyncratic state details 
allocate(e2Grid(ni), e3Grid(ni))
allocate(iTrans(ni, ni), iStat(ni))
allocate(yMat(ne1, ne2, ne3), yBarMat(ne1, ne2))
allocate(e3Map(ni), e2Map(ni))

! consumption and sunny day weights for non-default actions 
allocate(WD0(naNeg,ni,ne1), J7(naNeg,ni,ne1), J13(naNeg,ni,ne1), JR(naNeg, ni, ne1))
allocate(CA0(na, na, ni, ne1), dA0(na, na, ni, ne1), apiMax0(na, ni, ne1), apiMin0(na, ni, ne1), MustDefault0(na, ni, ne1), CantSave0(na, ni, ne1))
allocate(CA7(naPos, naPos, ni, ne1), dA7(naPos, naPos, ni, ne1), apiMax7(naPos, ni, ne1))
allocate(CA13(naPos, naPos, ni, ne1, nl), dA13(naPos, naPos, ni, ne1, nl), apiMax13(naPos, ni, ne1, nl))

! consumption, recovery, and recovery rates for default actions 
allocate(CD7(naNeg, ne1, ne2, ne3), xi7(naNeg, ne1, ne2, ne3), xiTilde7(naNeg, ne1, ne2, ne3))
allocate(CD13(ne1, ne2, ne3), xi13(naNeg, ne1, ne2, ne3), xiTilde13(naNeg, ne1, ne2, ne3), lambda13(naNeg, ne1, ne2, ne3), lpiMat(naNeg, ne1, ne2, ne3), lpwMat(naNeg, ne1, ne2, ne3))
allocate(xiR(naNeg, ni, ne1), xiTildeR(naNeg, ni, ne1))
allocate(pA0_GoDQ(naNeg+1, naNeg))
allocate(xiGrid(nXi), xiProb(nXi))

! calibration
allocate(ScaledInternalParameters(ncal))
allocate(Model_TargetMoments(ntarg))

! spreads by default probability calculations 
allocate(spreads(naNeg, na, ni, ne1), rates(naNeg, na, ni, ne1))
allocate(Cond_DefProb1(naNeg, ni, ne1))
allocate(Cond_DefProb(naNeg, ni, ne1, nAhead))

allocate(ExpDefProb1(na, ni, ne1))
allocate(ExpDefProb(na, ni, ne1, nAhead))
allocate(pDPartition(nPartition), pDPartitionPlot(nPartition))
allocate(AvgBinSpread(nPartition), StDevBinSpread(nPartition), MassPartition(nPartition))

! SOLVE THE MODEL 
print*, ' '
print*, ' '
print*, ' '
print*, 'SINGLE SOLUTION OF THE MODEL'
print*, ' '
print*, ' '
print*, ' '
call mapfull(ScaledInternalParameters, InternalParameters_LB, InternalParameters_UB, ncal, InternalParameters)
call SolveModelForGivenParameters(CriterionValue, Model_TargetMoments, ScaledInternalParameters)


end program SS_JPEM_Main